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f  generalized  nonsymmetric  eigenvalue 
matrices,  it  is  often  desirable  to  work 
*M  jib  order  to  enhance  the  eigenvalue 


''Abstract.  When  using  an  iterative  method  for  solving 
problem  of  the  form  Fu  —  \Mu,  where  F  and  M  aje  rea 
with  the  shifted  and  inverted  operator  B  —  (K  -  <tM)~ 
separation  and  improve  efficiency.  Unfortunately,  the  shift  4  is  generally  complex  and  so  is  the 
matrix  B.  The  question  then  is  whether  it  is  possible  to  avoid  complex  arithmetic  while  preserving 
the  advantages  of  bandedness  of  the  pair  ( F,M ).  For  the  classical  problem  where  M  =  I  and  F 
is  banded,  complex  arithmetic  can  be  avoided  by  using  double  shifts,  i.e.,  by  working  with  the  real 
matrix  BB  whose  bandwidth  is  double  that  of  F.  This  satisfactory  solution  extends  to  the  case 
where  M  is  diagonal  as  well.  In  the  generalized  case  the  answer  to  the  above  question  is  negative, 
in  the  sense  that  complex  arithmetic  can  be  avoided  only  at  the  expense  of  loosing  the  advantage  of 
bandedness.  One  solution  is  to  factor  the  shifted  matrix  F  —  bM  in  complex  arithmetic  but  employ 
real  arithmetic  subsequently  in  the  iterative  procedure.  This  paper  examines  several  approaches 
and  discusses  their  respective  merits  under  different  circumstances. 
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1.  Introduction 

This  communication  is  concerned  with  the  eigenvalue  problem:  solve 

(F  —  \M)z  =  0,  (1.1) 

for  A  £  C  and  2  €  C^,  when  F,M  are  real  N  x  N  matrices.  Throughout  the  discussion,  M  will  be 
symmetric  and  positive  definite.  In  a  number  of  applications,  M  =  /,  the  identity  matrix.  There 
will  be  no  restriction  on  F  except  that  it  have  complex  eigenvalues.  More  precisely,  we  suppose 
that  only  a  few  of  the  eigenvalues  A  are  wanted,  namely  those  in  the  vicinity  of  a  given  complex 
number  a.  What  makes  this  task  challenging  is 

(I)  the  desire  to  keep  computation  in  R  rather  than  in  C; 

(II)  the  desire  to  exploit  any  narrow  band  structure  enjoyed  by  F  and  M. 

These  two  desires  can  be  in  conflict.  We  shall  assume  for  convenience  that  F  and  M  have  the 
same  bandwidth  20  +  1;  the  (i,j)  elements  vanish  whenever  |t  -  j\>0.  Moreover,  we  shall  assume 
that  the  band  is  narrow,  i.e.,  that  0  <<  n.  Note  that  the  goal  (II)  can  be  generalized  into  that  of 
exploiting  any  particular  sparse  structure,  not  just  bandedness. 

Every  reasonable  approach  known  to  us  requires  an  iterative  process  at  each  stage  of  which  a 
system  of  equations  must  be  solved.  The  simplest  of  these  is 

(F  -  cM)y  =  Mb 

where  b  is  given  and  y  is  to  be  computed.  Our  problem  reduces  to  an  attempt  to  reconcile  the  two 
aims  (I)  and  (II),  when  solving  the  above  system,  or  rival  ones  similar  to  it. 

In  the  body  of  this  paper  we  present  all  the  alternatives  that  have  occured  to  us  and  analyze 
them.  In  particular,  we  show  a  surprising  connection  between  two  of  them.  Unfortunately,  our 
analysis  leads  to  no  “best”  method,  but  we  give  operation  counts  and  storage  costs  for  the  better 
techniques. 

Before  proceeding  with  the  algebra,  we  should  say  something  about  complex  arithmetic.  We 
can  imagine  an  arithmetic  engine  that  would  employ  4  real  arithmetic  processors  in  parallel  to 
compute  the  product  of  two  complex  numbers  in  almost  the  same  time  as  required  for  a  real 
multiplication.  We  know  of  no  such  computer  at  present.  In  some  systems  we  know  of,  (VAX  780), 
the  ratio  complexrreal  arithmetic  is  nearly  four  but,  in  others,  the  cost  of  accessing  the  arguments 
has  become  sufficiently  large  to  reduce  the  ratio  to  nearly  two.  The  storage  penalty  remains  at  2:1. 

2.  Inverse  Iteration 

Throughout  the  paper  we  write  the  shift  a  as 

<x  —  p  +  i0, 

with  0>O  and  i2  =  —  1.  In  our  context,  inverse  iteration  is  defined  as  follows: 

1.  Choose  satisfying  ||z^||  =  1. 

2.  For  k  =  1,2...,.  until  convergence  do 

2.1.  Solve 

[F -oM)yw  =  MxW,  (2.1) 

2.2.  Normalize:  =  y(fc)/£(*),  where  is  the  component  of  yW  of  largest  modulus. 

2.3.  Check  for  convergence. 
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In  the  generic  case  the  sequence  {x^}  converges  to  z  where  Fz  =  \Mz  and  A  is  the  eigenvalue 
closest  to  a.  One  can  approximate  A  by  tr+x^(j)/y^{j)  where  y^(j)  is  an  above  average  element 
of  A  better  but  more  expensive  approximation  is  the  Rayleigh  Quotient, 

(Fx<*>,s 
(Mx(k),x(k)) 

Of  course  we  may  seek  several  eigenvalues  close  to  cr,  not  just  one.  Consequently,  more  elaborate 
iterations  are  needed.  Examples  are  Arnoldi  [1,  6],  Lanczos  [4,  2],  or  simultaneous  iteration  [3] 
(also  known  as  subspace  iteration).  The  differences  between  these  methods  are  not  important  here 
because  they  may  all  be  used  with  the  same  operator,  namely 


B  =  (F  -  oM)~l M.  (2.2) 

Note  that  the  sequence  { x (*)}  may  be  thought  of  as  generated  by  multiplying  each  term  by  B 
and  then  mormalizing  in  order  to  get  the  next  one.  The  matrix  B  is  not  formed  explicitly.  The 
dominant  part  of  each  step  in  any  of  the  iterative  methods  is  the  solution  step  2.1  of  the  algorithm. 
One  way  to  carry  this  out  is 

Method  1 

to  compute  the  triangular  factorization, 

F-aM=  LU , 

once  and  for  all  (in  complex  arithmetic).  Then  system  (2.1)  is  solved  by  the  two  triangular 
solves  (in  complex  arithmetic): 

UyW  =  w{k). 


Recalling  the  two  goals  in  the  introdution,  we  see  that  by  abandoning  (I),  real  arithmetic,  we 
can  exploit  (II),  band  structure.  The  costs  of  Method  1  are  as  follows. 

•  Arithmetic.  Factorization  :  02N  complex  multiplications.  Forward  and  backward  solutions: 
2 (0  +  1  )N  complex  multiplications.  Normalization:  jV  comparisons  and  N  real  multiplications 
in  step  2.2. 

•  Storage.  2(20+1  )N  real  locations  for  F  and  M;  (2/b  +  l )N  complex  locations  for  F~aM  =  LU; 
plus  two  complex  vectors  for  storing  x and  y^. 

Now  consider  the  implementation  of  step  2.1  in  real  arithmetic.  We  write  y  =  yr  +  i  y, ,  for 
any  vector  y  6  CN .  In  the  standard  way  we  equate  real  and  imaginary  parts  to  get 

(F  -  pM)yr  +  dMyi  =  Mxr 

-0Myr  +  (M  -  pM)yi  =  Mx,- 

( F  -  pM  6M  \  (yr\  (M  O  Wxr\ 

^  -0M  F-pMj  {yj  VO  M )  \  x,  /  ’ 
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(2.3) 

(2.4) 

(2.5) 
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or,  in  matrix  form 


This  system  of  order  2 JV  has  lost  the  band  structure  of  (2.1).  It  achieves  (I)  at  the  expense  of 
(II).  For  future  reference  it  is  important  to  observe  that  the  iteration  associated  with  (2.5)  attempts 
to  compute  the  eigenpairs  of  the  (2iV)  x  (2ZV)  real  matrix 

(  M  0\~l  (F~pM  9M  \_(M~lF-pI  91  \ 

\0  M)  \  ~dM  F-pMj~\  -91  M~lF  -  pi )  ' 

The  expression  of  the  inverse  of  the  above  matrix  is  of  great  help  when  establishing  relationships 
between  the  various  approaches  taken  in  later  sections.  Letting 

A  =  M~lF 

we  have 

(A- pi  91  \X  _(X~'  -9Y~'\  . 

\  -91  A- pi)  ~  [OY-1  X-1  )’  l2,6' 

in  which, 

X=(A-pI)  +  92(A-pI)~\  Y  =  X(A  -  pi)  =  (A  -  pi)2  +  02I  (2.7) 

In  particular  it  can  be  seen  from  above  that  yr,  yi  can  be  obtained  by  solving  for  yT  first,  by 

[(F  -  PM)  +  92M(F  -  pM)~lM ]  yr  =  Mxr  -  9M(F  -  pM)~lMxi  (2.8) 

and  then  getting  y;  by  substitution  in  the  equation  (2.3), 

9Myi  =  Mxr  -  [F  -  pM)yr, 

which  gives 

Vi  =  \[xr~  {M~l  F  -  pl)yr] .  (2.9) 

Note  that  one  can  also  compute  y,  first  from  an  equation  similar  to  (2.8)  and  then  substitute  in 
(2.4)  to  get  yr. 

When  M  =  I  a  simplification  is  possible,  by  multiplying  both  sides  of  (2.8)  by  (F  —  pi): 
Method  2  (For  the  case  M  =  I) 

To  solve  the  system  (2.1)in  inverse  iteration  algorithm,  compute  the  real  part  yr  of  yW  by 

[{F  -  pi)2  +  92]  yr  =  (F  -  pl)xr  -  9Xi ,  (2.10) 

and  its  imaginary  part  y,-  by  (2.9). 

The  bandwidth  has  been  doubled  but  not  ruined.  The  matrix  F2  may  be  stored  once  and  for 
all,  allowing  for  changes  in  cr.  The  costs  of  this  method  are  as  follows. 

•  Arithmetic.  Factorization  of  (F  -  pi)2  +  92I  (done  once)  :  4/?2iV  real  multiplications.  Forward 
and  backward  solutions:  8 (/?  +  1  )N  real  multiplications.  Normalization:  2 N  comparisons  and 
2N  real  multiplications  in  step  2.2. 

•  Storage.  (40  +  l)iV  real  locations  (F  -  pi)2  +  92I,  and  (40  +  1) JV  real  locations  for  its  LU 
factorization.  Four  real  vectors  for  storing  x W  and  yW. 


3.  The  general  case. 

Consider  again  (2.8).  If  M  is  diagonal  then  Method  2  extends  readily  and  we  will  not  con¬ 
sider  this  further.  We  now  take  up  the  general  case  when  M  has  the  same  band  structure  as  F. 
Premultiply  (2.8)  by  ( F  -  pM)M~l  to  get  the  analogue  of  (2.10): 

[(F  -  pM)M~1  (F  -  pM )  +  92M]  yr=(F-  pM)xr  -  6Mx,  (3.1) 

We  will  define 

G  s  (F  -  pM)M~1{F  -  pM )  +  92M  =  FM~1F  -  2 pF  +  \a\2M.  (3.2) 

Unless  M  is  diagonal,  the  presence  of  A/-1  ruins  the  bandedness  of  the  matrix  G  which  will  be  full 
in  general.  Note  that 

M~lG  =  Y  (3.3) 

where  Y  is  defined  in  (2.7). 

Nevertheless,  band  structure  may  still  be  exploited,  especially  in  a  number  of  applications 
where  F  and  M  are  generated  by  the  finite  element  method.  In  those  situations  it  is  common  to 
replace  the  consistant  mass  matrix  M  by  a  diagonal  lumped  mass  marix  D.  The  diagonal  elements 
of  D  are  the  elements  of  the  vector  Me  where  e  =  (1,1,...,  1)T.  The  real  matrix 

FD~lF  ~1pF  +  \<j\2M 

has  twice  the  bandwith  of  F  and  M  and  may  be  factored  efficiently  into  LU.  This  matrix  is  used  as 
a  preconditioner  for  the  proper  matrix.  The  inner  iteration  should  converge  in  a  very  small  number 
of  iterations.  The  following  iteration  on  the  residuals  is  the  simplest  technique. 

Method  3 

Set  r  *—  (F  -  pM)br  -  0A/6,- 
Until  convergence  do: 

(i)  Solve  LUd  =  r 

(ii)  compute  r  *—  r  —  Gd,  yr  —  yr  +  d. 

(iii)  If  || r||  too  large  then  repeat. 

(iv)  else  get  y,-  by  (2.9)  and  return. 

We  omit  to  give  the  details  on  the  costs  of  the  above  method,  because  the  process  is  iterative 
and  is  not  comparable  to  previous  techniques. 

4.  The  double  shift  approach 

A  problem  similar  to  ours,  but  in  the  context  of  the  QR  algorithm,  was  solved  by  J.G.F. 
Francis  in  1961/62.  If  A  €  RjV,lV  and  a  €  C  then  Y  =  [A- <rI)(A-  9l)  =  [(A  —  pi)2  +  92I]  €  R^’^. 
This  matrix,  which  is  real,  is  a  quadratic  polynomial  in  A  and  shares  .4's  eigenvectors.  When  A 
is  replaced  by  M~lF  then  this  matrix  coincides  with  the  matrix  T  defined  by  (2.7).  By  (3.3),  the 
eigenvalues  of  Y  are  those  of  the  generalized  eigenvalue  problem  Gz  =  vz  or: 

[FM~lF  -2pF  +  \o\2M]z  =  uMz  (4.1) 

where  v  =  A2  -  2pA  +  |cr(2.  Unless  M  is  diagonal,  or  block  diagonal,  this  matrix  is  real  but  full. 
Even  the  actual  computation  of  this  matrix  may  not  be  practically  feasible. 
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In  the  present  context,  inverse  iteration  is  defined  as  in  section  2  except  that  the  system  (2.1) 
is  replaced  by 

Gyw  =  Mx(k\  (4.2) 

We  can  solve  the  above  system  iteratively  as  indicated  in  the  previous  section.  However,  there 
are  alternative  approaches  that  fully  exploit  band  structure,  provided  we  relax  our  constraint  on 
working  entirely  in  real  arithmetic.  Let  LU  be  the  (complex)  factorization  of  ( F  -  <j M).  The 
matrices  L  anf  U  will  inherit  the  band  structure  of  F  and  M.  Now  to  solve  Gy  =  M x  where 
y  €  R^,  x  E  R^,  note  that 

Gy=LUM~lWy  =  Mx. 

An  algorithm  for  computing  y  is 

Method  4 

(i)  Solve  La  =  Mx  for  a  e.  C'v 

(ii)  Solve  Ub  =  Ma  for  ab  E  CN 

(iii)  Form  c  =  Mb 

(iv)  Solve  Ld  =  c  for  d  E  CN 

(v)  Solve  Ue  =  d  for  e  E  CN 

(vi)  Sety=  Re  (e). 

The  complex  arithmetic  is  hidden  in  the  above  subroutine  that  maps  x  into  y.  The  iteratiom  that  is 
used  to  compute  one  and  two  dimensional  eigenspaces  of  Y~l  can  confine  itself  to  real  arithmetic. 
We  shall  have  more  to  say  about  the  matrix  Y~l  in  the  next  section.  Now  we  resume  the  quest  for 
an  operator  that  requires  no  complex  arithmetic  and  yet  takes  advantage  of  narrow  bandwidth. 

5.  Real  and  Imaginary  Part  approaches 

Inverse  iteration  with  shift  a  is  equivalent  to  direct  iteration  (i.e.,  the  power  method)  using 
the  operator  ( F  -  crM)~l  on  C^.  To  obtain  related  operators  on  R^  we  can  take  the  real  and 
imaginary  parts 


B+  =  l  [(^  -  tfM)-1  +  {F  -  dM)~l]  M  ~  Re  [(F  -  oM)~'M] , 

it 

B-  =  ^  [(F  -  <rM)~ 1  -  (F  -  ffAf)-1]  M  =  Im  [(F  -  . 

If  Fz  =  \Mz  then 


B+z 


1 

2 


'  1  1 
A  —  a  A  -  9 


B~z 


1_ 

2i 


A  -  a 


(5.1) 

(5.2) 


(5.3), 


defining  and  the  eigenvalues  of  B+  and  B-  associated  with  the  eigenvector  z  of  A.  It  is 
readily  verified  that  as  A  — +  <r, 


H±  * 


1 


2(A  -  <t)  ' 


Thus  B+  and  B _  give  the  same  enhancement  to  eigenvalues  close  to  zero.  In  contrast,  as  A  — ►  oc, 
B-  dampens  the  eigenvalues  more  strongly  than  does  B+  since, 


A  -  p  _  9 

^  ~  (, A  — <r)(A-*)’  “  (A  -  cr)(A  -9)' 


(5.4) 
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The  question  now  is  whether  it  is  possible  to  reconcile  the  two  aims  (I)  and  (11)  set  out  in  the 
introduction,  when  computing  B+v  and  B-v  for  any  v  €  R^.  Reference  back  to  (2.8)  shows  that 

B+v  =  [(F  -  pM )  4-  92M{F  -  pM)~xM]~l  v  =  X~xv. 

If  real  arithmetic  is  mandatory  (aim  (I))  then  the  presence  of  the  full  matrix  ( F  -  pM)~x  in  X 
precludes  the  exploitation  of  bandedness  (aim  (II))  in  the  triangular  factorization  of  X.  This  leaves 
two  possibilities: 

1.  Solve  -Yu  =  v  for  u  iteratively,  in  real  arithmetic,  exploiting  band  structure  as  described  at  the 
end  of  Section  3. 

2.  Ignore  the  structure  of  B+  and  evaluate  B+v  by  solving  (F-crM)u  =  Mv  in  complex  arithmetic 
and  then  returning  the  real  part  (respectively  the  imaginary  part  for  a  method  using  B_  )  This 
is  Method  5.  • 

Method  5 

1.  Solve  (F  —  <jM)w  =  Mv  (  complex  arithmetic). 

2.  Set  B+v  —  Re(w)  (respecively  B-v  =  Im(w)). 

The  cost  of  the  above  method  is  as  follows. 

•  Arithmetic.  Factorization  (done  only  once)  :  (32N  complex  multiplications,  Forward  and  back¬ 
ward  solutions:  2(/3+l)iV  complex  multiplications.  Normalization:  N  comparisons  and  N  real 
multiplications  in  step  2.2. 

•  Storage.  2(2(3 +l)N  real  locations  for  F  and  M\  (2/9+1) .V  complex  locations  for  F-aM  =  LU\ 
plus  two  complex  vectors  for  storing  and  y(kV 

Method  5  is  a  compromise.  What  must  be  emphasized  here  is  that  from  the  point  of  view  of  the 
iterative  methods,  such  as  Arnoldi,  Lanczos,  or  subspace  iteration,  that  will  be  making  use  of  B+ 
there  is  no  compromise.  Goal  (I)  is  realized.  These  iterations  will  use  real  arithmetic  exclusively. 
Goal  (II)  is  achieved  by  using  complex  arithmetic  in  the  lower  level  subroutine  that  evaluates  B+v. 

Note  that  the  cost  of  Method  5  is  lower  than  that  of  Method  4  of  the  previous  section.  In  fact 
the  extra  work  in  Method  4  brings  no  further  benefit  in  the  light  of  the  following  surprising  result. 
Recall  that  G  =  (F  -  aM)M~x(F  -  aM)  +  92M. 

Theorem  5.1.  The  matrices  B-,G  and  M  are  related  by 

B _  =  9G~lM.  (5.5) 


Proof.  We  have,  by  definition, 

B-  =  j.  [(F  -  aM)-1  -  (F  -  oM)~x]M 

=  i(F  -  [(F  -  9M)  -  (F  -  crM)\  (F  -  9M)~lM 

=  (F  -  crM)-x0M(F  ~  9M)~xM 
=  9G~XM 
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By  the  theorem  the  solution  y^  of  (4.2)  is  identical  with  B~x^k\  apart  from  the  multiplicative 
scalar  9.  Note  also  that  the  right-hand  side  of  (5.5)  is  nothing  but  9Y~l ,  i.e.,  the  block  in  position 
(2,1)  of  the  matrix  in  (2.6). 


6.  Numerical  experiments 

All  numerical  tests  have  been  performed  on  a  Vax-785  using  double  precision,  i.e.,  the  unit 
roundoff  is  2-56  as  1.3877  x  10~17.  Our  test  example,  taken  from  [5],  models  concentration  waves 
in  reaction  and  transport  interaction  of  some  chemical  solutions  in  a  tubular  reactor.  The  concen¬ 
trations  x(t, z),y(r,z)  of  two  reacting  and  diffusing  components,  where  0  <  z  <  1,  represents  a 
coordinate  along  the  tube,  and  where  r  is  the  time,  are  modeled  by  the  system  [5]: 


dx  _  Dz  drx 
dr  L 2  dz 2 

dy_  _  Dv  d2y 
dr  L2  dz 2 

with  the  initial  condition 


/(*,  y). 

(6.1) 

g{x,y), 

(6.2) 

x(0,z)  =  x0{z),  y{0,z)  =  y0(z),  'i  z  e  [0,1], 

and  the  Dirichlet  boundary  conditions: 

z(0,r)  =  x(l,r)  =  x *  ,  y(0,r)  =  y(l,r)  =  y*. 

We  consider  in  particular  the  so-called  Brusselator  wave  model  [5]  in  which 

f[x,y)  -  A  -  (B  +  l)x  +  x2y  ,  g(x,y)  =  Bx  -  x2y.  (6.3) 


Then,  the  above  system  admits  the  trivial  stationary  solution  a:*  =  A,  y*  =  B/A. 

In  this  problem  one  is  primarily  interested  in  the  existence  of  stable  periodic  solutions  to  the 
system  as  the  bifurcation  parameter  L  varies.  This  occurs  when  the  eigenvalues  of  largest  real  parts 
of  the  Jacobian  of  the  right  hand  side  of  (6.1)  -  (6.2),  evaluated  at  the  steady  state  solution,  is 
purely  imaginary.  For  the  purpose  of  verifying  this  fact  numerically,  one  first  needs  to  discretize 
the  equations  with  respect  to  the  variable  z  and  compute  the  eigenvalues  with  largest  real  parts  of 
the  resulting  discrete  Jacobian. 

The  exact  eigenvalues  are  known  and  this  problem  is  analytically  solvable.  The  article  [5] 
considers  the  following  set  of  parameters 

Dz  =  0.008,  Dv  =  \dx  =  0.004,  A  =  2,  B  =  5.45. 

For  small  L  the  Jacobian  has  only  eigenvalues  with  negative  real  parts.  At  L  as  0.51302  a  purely 
imaginary  eigenvalue  appears. 

We  discretize  the  interval  [0, 1]  using  n  interior  points,  and  define  the  mesh  size  h  =  l/(n+  1). 
The  discrete  vector  is  of  the  form  where  x  and  y  are  n-dimensional  vectors.  We  denote  by 
fh  and  gh  the  corresponding  discretized  functions  /  and  g ,  the  Jacobian  is  a  2  x  2  block  matrix  in 
which  the  diagonal  blocks  (1, 1)  and  (2,2)  are  the  matrices 
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Tridiag{l,  -2, 1}  + 


dhjx,  y) 
dx 
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and 

Tridiag{l,~ 2, 1}  + 

respectively,  while  the  blocks  (1,2)  and  (2,1)  are 


dgh{x,y) 

dry 


dfh{x,y) 

dy 


and 


d9h{*,y) 

dx 


respectively.  Note  that  since  the  two  functions  /  and  g  do  not  depend  on  the  variable  z,  the 
Jacobians  of  either  fh  or  gh  with  respect  to  either  x  or  y  are  scaled  identity  matrices.  We  denote 
by  A  the  resulting  2n  x  2n  Jacobian  matrix.  In  the  following  tests  we  took  n  =  100,  which  yields 
a  matrix  A  of  size  200.  We  point  out  that  the  exact  eigenvalues  of  A  are  readily  computable,  since 
there  exists  a  quadratic  relation  between  the  eigenvalues  of  the  matrix  A  and  those  of  the  classical 
difference  matrix  T ridiag{l,  -2, 1}.  In  fact  part  of  the  spectrum  (the  32  rightmost  eigenvalues)  of 
the  matrix  A  is  shown  in  Figure  4.  We  have  not  shown  the  rest  of  the  spectrum  of  A  consisting 
of  168  real  eigenvalues  that  are  almost  uniformely  distributed  in  the  interval  [-1,235.5,-51.912]. 
The  rightmost  eigenvalues,  determined  with  maximum  accuracy,  i.e.,  approximately  16  digits  are 

A1i2  =  1.8199876787305946  x  10-5  ±  2.139497522076329 


-\s  is  observed  the  real  part  is  close  to  zero,  which  verifies  the  theory,  within  discretization  errors. 

The  purpose  of  these  experiments  is  to  compare  the  performances  of  the  methods  using  the 
three  approaches  B  =  (A  — <r/)-1  ,  B+  =  Re(B)  and  B-  =  Im(B),  all  in  conjunction  with  Arnoldi’s 
method.  We  have  plotted  the  convergence  history  for  the  three  methods  for  three  choices  of  the 
shift  cr,  namely  a  =  0.1  +  2. 1 1,  a  —  0.0  +  2.5* ,  and  a  —  0.5  +  2. It.  The  plots  in  Figures  1,  2  and  3 
show  the  relative  errors 


|  Ar  | 

versus  the  number  of  Arnoldi  steps.  As  is  observed  the  performances  of  the  two  different  approaches 
are  not  constant. 

In  Figures  5  and  6  we  show  the  spectra  of  the  corresponding  matrices  J3+  and  B_  for  the  last 
two  cases,  i.e.,  for  a  =  0.5  +  2.1*  and  for  cj  =  2.5i.  In  each  case  we  have  circled  the  eigenvalue  of 
largest  modulus.  Notice  the  very  good  separation  properties  of  the  dominant  eigenvalue  despite  a 
relatively  distant  shift.  Also  observe  the  concentration  around  the  origin  of  the  transformed  large 
eigenvalues  of  A.  The  reader  should  note  that  the  scales  are  different.  For  example  in  the  top 
graph  in  Figure  6  the  x-axis  has  a  total  length  of  0.08,  which  means  that  the  spectrum  is  almost 
purely  imaginary  in  this  case.  The  spectra  of  B+  and  B_  bear  no  particular  resemblance  and  it  is 
hard  to  predict  from  looking  at  the  pictures  only  which  method  will  converge  faster. 

It  is  also  instructive  to  compare  the  two  mappings  ^+(A)  and  /*_(A)  as  defined  by  (5.4).  .As 
an  experiment  we  plotted  the  images  /u+(A)  and  /*-(A)  of  several  circles  of  small  radii,  centered  at 
the  shift  a.  The  goal  is  to  compare  the  two  mappings  for  a  similar  situation  where  the  eigenvalues 
of  the  original  matrix  A  are  distributed  in  circles  aiound  the  shift.  When  a  =  0.5  +  2.1*  and  the 
radii  were  0.1, 0.2, ...  ,0.5  respectively,  the  two  resulting  plots  looked  very  much  like  five  concentric 
circles  and  were  almost  indistinguishable  for  the  cases  B+  and  B-.  For  this  reason  we  omit  to  show 
the  resulting  figures.  Changing  a  to  0.5  and  taking  the  same  radii  as  above  produced  the  graphs 
in  Figure  7. 

Notice  again  that  the  outmost  curves,  those  corresponding  to  the  dominant  eigenvalues  of  B+ 
and  B_,  are  slight  perturbations  of  circles.  Although  not  apparent  at  first  glance,  the  outmost 
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curve  for  B+  is  almost  superposable  with  that  of  B-  provided  we  shifted  the  whole  plot  of  B- 
in  the  south-west  direction  by  about  one  unit  of  the  graph.  Notice  also  that  the  B-  graph  is 
symmetric  about  the  real  axis  as  is  expected  from  the  definition  of  /z_.  Similarly,  the  n+  curves 
can  be  seen  to  be  symmetric  with  respect  to  the  imaginary  axis,  as  is  verified  in  the  plots. 

In  the  following  discussion  we  assume  that  there  is  one  actual  eigenvalue  of  A  per  circle, 
i.e,  there  is  only  one  eigenvalue  of  B+  or  B-  per  curve  represented  in  the  two  plots  of  Figure  7 
repectively.  .Assume  at  first  that  the  two  dominant  eigenvalues  for  B+  are  located  on  the  imaginary 
axis  in  the  lower  half  plane.  These  are  roughly  (i+, i  ^  -5.45 i  and  /r+i2  %  -2.9 i  which  means  that 
the  the  convergence  ratio  in  inverse  iteration  would  be  ^+,i\  ^  0.532.  It  is  found  that  the 

corresponding  eigenvalues  of  the  £?_  matrix  are  the  two  dominant  eigenvalues  -4.54  and  -2.08. 
which  leads  to  the  convergence  ratio  of  0.402,  much  better  than  that  of  the  B+  approach.  Assume 
on  the  rther  hand  that  both  the  dominant  and  the  subdominant  eigenvalues  of  B-  are  located  on 
the  real  axis  on  the  right  half  plane:  /z_it  «  5.554  and  /i_, 2  ^  3.12.  Then  the  associated  convergence 
ratio  for  inverse  iteration  with  B_  becomes  ,1 1  %  0.562.  The  corresponding  eigenvalues 

of  B+  are  found  to  be  approximately  ji+ii  4.44  and  ft*  1.875  which  gives  the  convergence 
ratio  0.44  for  inverse  iteration  with  £+,  a  much  better  ratio  than  that  of  the  £_  approach.  Thus, 
the  previous  situation  has  been  completely  reversed.  What  is  interesting  is  that  this  has  occured 
in  spite  of  keeping  the  distances  of  the  two  eigenvalues  of  A  closest  to  the  shift  the  same  in  both 
situations.  In  other  words  it  is  not  only  the  distance  of  these  eigenvalues  that  matters  for  the  speed 
of  inverse  iteration,  but  also  their  relative  location  around  the  shift.  This  tells  us  that  in  practice 
it  will  be  vain  to  try  determining  a-priori  which  of  the  two  approaches  is  to  be  favored. 

7.  Summary  and  conclusion 

We  have  examined  several  ways  of  implementing  shift  and  invert  techniques  for  the  eigenvalue 
problem  F:  =  XMz ,  in  the  common  situation  where  the  shift  is  complex  while  F  and  M  are  real  and 
banded  matrices.  If  M  is  diagonal  there  are  several  possible  variants  which  will  perform  equally 
well.  On  the  other  hand,  when  M  is  arbitrary,  then  any  attempt  to  avoid  complex  arithmetic 
completely  would  result  in  problems  with  full  matrices.  Then  the  advantages  that  might  be  gained 
from  any  exploitable  structure  of  F  and  M,  such  as  bandedness  sparsity  and  so  on,  would  be  lost. 
One  alternative  proposed  for  this  situation  is  to  use  a  real  operator,  whose  eigenvalues  offer  the  same 
separation  enhancement  as  those  of  the  ideal  operator  B  =  (F  —  <rM)~l M .  Two  such  operators  are 
the  real  part  B+  and  the  imaginary  part  B _  of  the  operator  B.  Thus,  in  a  typical  iterative  method, 
for  example  Arnoldi’s  method,  the  factorization  of  the  operator  B  and  the  forward  and  backward 
solutions  needed  when  applying  B  to  a  vector,  are  still  performed  in  complex  arithmetic,  but  the 
iterative  method  itself,  e.g.  Arnold!,  would  be  realized  entirely  in  real  arithmetic.  The  first  dvantage 
of  using  either  of  these  approaches  over  that  of  using  B,  is  economical:  we  save  computational  time 
and  storage,  in  the  Arnoldi  part  of  the  method.  The  second  is  purely  mathematical:  we  have 
replaced  an  eigenvalue  problem  with  a  real  operator  with  one  having  the  same  property. 

Although  it  is  practically  infeasible  to  determine  which  of  the  two  approaches  B+  or  B-  is 
best  in  general,  the  numerical  experiments  suggest  that  they  both  perform  nearly  as  well  as  that 
of  using  the  operator  B. 
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